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ABSTRACT 

Context. Forthcoming experiments will enable us to determine tomographic shear spectra at a high precision level. Most predictions 
about them have until now been biased on algorithms yielding the expected linear and non-linear spectrum of density fluctuations. 
Even when simulations have been used, so-called Halofit (Smith et al 2003) predictions on fairly large scales have been needed. 
Aims. We wish to go beyond this limitation. 

Methods. We perform N-body and hydrodynamical simulations within a sufficiently large cosmological volume to allow a direct 
connection between simulations and linear spectra. While covering large length-scales, the simulation resolution is good enough to 
allow us to explore the high-^ harmonics of the cosmic shear (up to I ~ 50000), well into the domain where baryon physics becomes 
important. We then compare shear spectra in the absence and in presence of various kinds of baryon physics, such as radiative cooling, 
star formation, and supernova feedback in the form of galactic winds. 

Results. We distinguish several typical properties of matter fluctuation spectra in the different simulations and test their impact on 
shear spectra. 

Conclusions. We compare our outputs with those obtainable using approximate expressions for non-linear spectra, and identify 
substantial discrepancies even between our results and those of purely N-body results. Our simulations and the treatment of their 
outputs however enable us, for the first time, to obtain shear results taht are fully independent of any approximate expression, also 
in the high-f range, where we need to incorporate a non-linear power spectrum of density perturbations, and the effects of baryon 
physics. This will allow us to fully exploit the cosmological information contained in future high-sensitivity cosmic shear surveys, 
exploring the physics of cosmic shears via weak lensing measurements. 

Key words. Cosmology: theory, Dark Matter and Dark Energy, gravitational lensing: weak, surveys 

1. Introduction (DEFT: Albrecht et al. 2006), while cosmic shear measures 

_ , ^„ s . , , , _ ,. , , have already been performed by using either large area surveys 

Dark energy (DE) is the most remarkable finding but largest carried Qut ^ nd _ based te lescopes (e.g., Hoekstra et al. 

outstanding uncertainty in cosmology today. Cosmic microwave 2QQ6 Fu et ^ 2m) Qr namjwer ^ g s wkh M h H 

M • background (CMB) spectra firmly constrain its contribution to from HST Ma ^ ^ Schrabback et al. 

the cosmic mass density budget, but its equation of state, w(z) 2010) observational campaigns to ultimately perform the sys- 

can be constrained only by measures of the mass density field tematic { of t0 hic cosmic shear have been 

p(x, z), at low z It is therefore important that the analysis of to- d for fumre missionSi such as the recentl approved E SA 

orographic weak lensing spectra, which are directly sensitive to Eudid satelHte (Laureijs et aJ 2(m) 
the whole matter distribution, can be suitably translated into in- 
formation on the density fluctuation spectrum 

P(k 7) = (Wk zll 2, > (11 ^ n essent ial point is that, according to Huterer & Takada 

(2005), spectral predictions suitable to exploit lensing data need 

and its redshift dependence. Here 6(k, z) is the Fourier transform a precision 0(1 %) in P(k, z). In turn this means that cosmolog- 

of the matter fluctuation field e(x, z) - p(x, z)/p- 1 . ical parameter variations causing a shift 0( 1 %) in P(k, z) could 

This paper is therefore devoted to examining the relation be- become appreciable. More recently, Hearin et al. (201 1) recon- 

tween P{k, z) and the tomographic weak lensing spectra Pij(€), sidered the whole question, showing that the precision require- 

which are defined below. This possibility has been debated by ments depend on ( max (the highest multipole included in the 

various authors and discussed by the Dark Energy task force lensing survey), the precision achieved in photometric redshift 
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(photo-z, hereafter) measurements, and the simultaneous use of 
other information, such as degree scale clustering data. 

The relation between P(k,z) and Pij{€) can be studied by 
designing simulations, and it is now clear that any attempt to 
measure the DE state equation w(a) requires spectral predic- 
tions even more precise than 1 %. Accordingly, we develop sim- 
ulations here to enable us (i) to connect their spectra on large 
scales with linear predictions, and (ii) to provide predictions up 
to ( =s 20000 - 30000, a range that experiments will test but 
no simulation ever yet explored without simultaneously using 
spectral approximations, which unavoidably affect the normal- 
ization of spectral predictions, even at large t. Our objective (i) 
is achieved by using a box of 410/r'Mpc on a side. Objective 
(ii) requires a sufficient large dynamical range, appropiate tech- 
niques to calculate the high-A: spectra, and an accurate evaluation 
of the numerical noise, which is initially present on the scale of 
the grid where the initial conditions are defined. 

The cosmology that we use to test this set of techniques 
is a Lambda cold dark matter (ACDM) model, for which the 
DE state equation is w(z) = -1 ■ We recall that, by using 
a set of ACDM N-body simulations in boxes of lengths of 
~ 240 /T 1 Mpc, consisting of 256 3 particles, Smith et al. (2002, 
see also Jenkins et al. 1998) produced numerical expressions 
for P(k,z), based on the halo model, claiming a precision of 
0(±3 %). These expressions, dubbed Halofit (hereafter HF), 
have also been inserted into linear calculations and, although 
imprecise even for ACDM (Casarini et al. 2009, Hilbert et al. 
2009, Heitmann et al. 2010), have been tentatively used to work 
out shear spectra even for cosmologies with DE state equations 
w + -1. As might be expected, this extension outside the as- 
signed range of validity causes misleading predictions, namely 
when parameter estimates, performed using tomographic shear 
spectra, are correlated (e.g. Casarini et al. 2011a, Seo et al. 
2011). 

Besides these problems, going beyond HF expressions is 
a must when exploring scales where baryon physics becomes 
important. In this paper values of k > 2-3/zMpc 1 are ex- 
tensively treated, as they correspond to a spectral range where 
cosmic shear will be accurately measured by future surveys. 
Accordingly, the simulations used in this work treat baryon 
physics by including adiabatic cooling, star formation, and en- 
ergy feedback from supernovae (SNe). Comparing them with 
simulations run with purely gravitational forces, or in the pres- 
ence of simplified baryon physics, will enable us to follow the 
impact of different physical effects on both fluctuation and cos- 
mic shear spectra and gauge the level of approximation in shear 
spectra when using HF. Several authors have dealt with these 
problems in the past few years. Jing et al. (2006) were the first to 
estimate the impact of radiative gas cooling on fluctuation spec- 
tra measured from hydrodynamical simulations, finding shifts of 
(9(10%) in the 1-10 liMpc' 1 fc-range. Rudd et al. (2008), us- 
ing a box of 60 h~ l Mpc and (2x)256 3 particles, compared spec- 
tra obtained using pure gravitational dynamics, with those based 
on either non-radiative or radiative baryon physics. Casarini et 
al. (2011a, 2011b) made use of boxes of 256 and 64/T 1 Mpc 
on a side, with (2x)256 3 particles, to analyse matter fluctua- 
tion spectra and shear spectra in models with different DE state 
equations, and studied the power of shear experiments to dis- 
criminate between them. Fedeli et al. (201 1) investigated the ef- 
fect of baryon physics on the matter power spectra for standard 
and extended quintessence models using a box of 300 /z _1 Mpc 
and (2x)768 3 particles. Viel et al. (201 1) explored the effects of 
baryons of the characteristic cut-off in the initial power predicted 
by the warm dark matter model using a set of simulations of size 



6.25-100 /r'Mpc and (2x)512 3 particles. We note also the study 
of Guillet et al. (2010), who did not derive spectra but, almost 
equivalently, studied the impact of gas dynamics on the variance 
(and skewness) of the mass distribution in a box of 50/z~'Mpc 
with (2x)1024 3 particles. Finally, Van Daalen et al. (2011) per- 
formed hydrodynamical simulations in boxes of various size, up 
to 100 h^ 1 Mpc on a side, by using (2x)5 1 2 3 particles and includ- 
ing AGN feedback, and used them to evaluate fluctuation spec- 
tra. In turn, Semboloni et al. (2011) used such spectra to derive 
their shear spectra. 

Our simulations are run in a box of side L = 410 h' 1 
Mpc (corresponding to ki ~ 1.5 x 10~ 2 /zMpc~', where spec- 
tra are perfectly linear even today) and have a force resolution 
e = 7.5 hT x kpc (corresponding to k e ~ 8.4 x 10 2 /zMpc _1 ), so 
that they can be used to produce data with high k values. More 
precisely, spectra can be generated up to k ^ N(2n/L) with N - 
2 15 = 32768, while maintaining k (^ 5.0 x lO^Mpc^ 1 ) < k e . 
A technique to obtain exact spectra up to such high-fe, without 
the use of exceedingly large grids, is described in Section 3. The 
spectra we mostly use extends up to k ^ 250 /zMpc -1 . 

The plan of the paper is as follows. We present in Section 
2 the simulations on which our analysis is based. In Section 3, 
we first describe the method to compute the power spectrum of 
density fluctuation from simulations. After discussing the effect 
of numerical noise and the comparing N-body results with those 
of HF, we then show the effect of baryon physics on the power 
spectrum. In Section 4, we pass from the fluctuation spectrum to 
the cosmic shear spectrum. We discuss our results and draw our 
main conclusions in Section 5. 

2. Simulations 

Our simulations follow the development of large-scale structure 
within a periodic box of comoving size L - 410/z~'Mpc, using 
(2x)1024 3 particles, in a spatially flat Gaussian ACDM model 
with Q m = 0.24, Q. h = 4.13 x 10~ 2 , h = 0.73, and n, = 0.96, 
(density parameters of total matter and baryons, Hubble param- 
eter, primordial spectral index, respectively) and are normalized 
so that erg = 0.8 at z = . With this parameter choice, our simu- 
lated cosmological model is consistent with the WMAP-7 CMB 
results (Komatsu et al. 2010). 

For the hydrodynamical simulations, we then have two pop- 
ulations of particles, whose mass ratio is chosen to reproduce 
the cosmic baryon fraction, initially placed on two uniform 
grids, displaced by half a grid size. Initial displacements from 
the unperturbed grid position are been generated according to 
the Zeldovich approximation at the initial redshift of z,„ =41. 
The resulting masses of CDM and baryonic particles are m c ^ 
1.89 X 10 9 h- l M Q and m b =* 3.93 x 10 8 h~ l M e , respectively. 

For the sake of comparison, we also run a purely gravita- 
tional simulation. In this case, the particles are initially set on a 
single uniform grid and have a mass m c ^ 2.28 x 10 9 /z _1 M o . 
They are then again displaced according to the Zeldovich ap- 
proximation to the initial redshift z,„. 

Simulations are carried out using the TreePM smoothed- 
particle hydrodynamics (SPH) GADGET-3 code, an improved 
version of the GADGET-2 code (Springel 2005), which al- 
lows each processor to also be assigned disjoint segments of 
the Peano-Hilbert curve, on which the decomposition of the 
computational domain is based. This provides a substantial im- 
provement to the balance among the work-loads assigned to 
the processors, thus substantially increasing the code efficiency. 
Gravitational forces are computed using a Plummer-equivalent 
softening, which is fixed to epi = 7.5 h physical kpc from z — 
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to z — 2, and fixed in comoving units at higher redshift. The ini- 
tial conditions are generated using the package NgenIC, part of 
the GADGET distribution. 

The simulation conteining only dark matter (DM) is dubbed 
DMo. We then have two hydrodynamical simulations. One is 
a previously developed non-radiative simulation (gravitational 
heating hereafter GH) that uses 64 neighbours to compute the 
hydrodynamic forces, with the width of the B-spline smoothing 
kernel allowing us to reach a minimum value equal to half of the 
gravitational softening. A second hydrodynamical simulation is 
carried out by including the effect of cooling and star forma- 
tion (hereafter CSF). In this simulation, radiative cooling is com- 
puted for non-vanishing metallicity according to Sutherland & 
Dopita (1993), who also considered heating/cooling from a spa- 
tially uniform and evolving ultraviolet background. Gas parti- 
cles above a given threshold density are treated as a multi-phase 
medium, so as to provide a subresolution description of the inter- 
stellar medium according to the model described by Springel & 
Hemquist (2003). Within each multi-phase gas particle, a cold 
and a hot-phase coexist in pressure equilibrium, with the cold 
phase providing the reservoir of star formation. Conversion of 
collisional gas particles into collisionless star particles proceeds 
in a stochastic way, with gas particles spawning a maximum 
of two generations of star particles. The CSF simulation also 
includes a description of metal production from chemical en- 
richment by supernova (SN) SNII, SNIa, and asymptotic giant 
branch stars, as described by Tornatore et al. (2007). Stars of 
different mass, distributed according to a Salpeter initial mass 
function (IMF), release metals over the timescale determined by 
the corresponding mass-dependent life times. Kinetic feedback 
is implemented by mimicking galactic ejecta powered by SN ex- 
plosions. In these runs, galactic winds have a mass upload pro- 
portional to the local star-formation rate. We use v w = 500km/s 
for the wind velocity, which corresponds to an assumption of 
about unity efficiency for the conversion of energy released by 
SNII into kinetic energy for a Salpeter IMF. The same sim- 
ulations were employed to analyse the contribution of baryon 
physics to the halo mass function (Cui et al. 201 1). 

The feedback model included in the CSF simulation is 
known not to be fully adequate for regulating overcooling, es- 
pecially in large cluster-sized halos (e.g., Borgani et al. 2006). 
It is also known that including AGN feedback leads to a more 
efficient regulation of star formation inside galaxy clusters. In 
turn, over-cooling is certainly absent in non-radiative simula- 
tions. Accordingly, CSF and NR simulations represent somehow 
two extremes, and a fair description of physical feedback mech- 
anisms is expected to yield the results derived between them. 
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Fig. 1. Upper panel: power spectra for the DM-only simulation, 
obtained by using n 3 grids with effective n from 1024 to 32768, 
are shown at z - (solid curves) and z - 4.62 (dashed curves). 
The latter is the highest redshift at which non-linear spectral fea- 
tures are (barely) visible above the numerical noise. In contrast, 
at z - 0, numerical noise affects only the spectrum at the highest 
resolution. Lower panel: The ratios of P/(k) (/=0,.., 4) to Ps(k) 
are plotted. Values of this ratio at the two considered redshifts 
are almost indistinguishable and this plot explicitly shows that 
they do not depend on the extension of numerical noise and are 
self-similar when the resolution limit is shifted. 



the density field is assigned reaches ^ 12.5 h~ l kpc. This folding 
technique was first proposed by Jenkins et al (1996) and a wider 
description can be found in Smith et al. (2003). In Appendix 
A, we briefly discuss why this tecnique is so effective. The ca- 
pability of this technique to yield spectra down to wavelengths 
slightly above the gravitational softening scale is limited by the 
numerical noise of the grid used to set the initial conditions. The 
problem, however, is not caused by the reductions in the box size 
and would be identical if large n grids could be directly applied 
to the original box of side L. 



3. Fluctuation spectra 

Power spectra are computed at 20 redshift values obtainable 
from the expression 1 + z r — 10 r/,2 ° (r - 0, 1, ..., 19). Spectra are 
evaluated by using the algorithm PMpowerM included in the PM 
package (Klypin and Holtzman 1997), courtesy of A. Klypin. 
Through a CiC procedure, the algorithm assigns the density field 
on a uniform cartesian grid starting from the particle distribution. 
It then evaluates the spectrum by applying Fast Fourier Tras form 
(FFT) on a n 3 grid. 

Here we consider the effective n - 2-OV values, with / from 
to 5, i.e. n from 1024 to 1024 x 2 5 = 32768 . Such large n 
are obtainable by considering a N 3 grid in a box of side Ljl) ", 
where all simulation particles are inset, in points of coordinates 
*i,f - M - vL/2? (i =1,2, 3), with an integer v selected so that 
< Xij < L/2-t . With this procedure, the size of the grid where 



3.1. Dark matter simulation: spectral extension and 
numerical noise 

In Figure Q] (upper panel), we overlap the spectra P/(k) of the 
DMo simulation at zo - and z\ - 4.62, obtained for values 
of the folding parameter / ranging from to 5, i.e. by using 
an effective n 3 grid, with n ranging from 1024 to 2 5 x 1024. At 
low k, baryonic acoustic oscillations (BAOs) are visible, with the 
expected amplitude, that are slightly attenuated from z\ to zq. 

While the onset of non-linearity is evident in the spectrum 
at z = 0, the situation is quite different at the higher redshift. 
Here numerical noise affects k values even below ~ 10, almost 
completely precluding the analysis of nonlinearity, displaying its 
main effects above such scale, at z\ ■ 

We note however that the shapes of the smearing of the Pj 
spectra (/ = ... 4), with respect to P$ , are unaffected by numeri- 
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Fig. 2. Smoothed (see text) logarithmic derivatives of P(k) at var- 
ious z values (here / = 5). 



cal noise. In the lower panel the cuts for zo and z\ are indeed both 
plotted, but they can hardly be distinguished. In the lower panel 
of Figure [T] one can also appreciate that the smearing shape is 
similar for any / value, just logarithmically displaced either up- 
or downwards. 

These results highlight that two criteria that need to be con- 
sidered when making an optimal choice of the number / of fold- 
ings, which is large enough to allow us to make use of spectra 
at large k values. At the same time, one should consider the on- 
set of numerical noise. In general, it is futile to attempt to derive 
spectra in regions where the numerical noise dominates. 

We then consider, for instance, the case n = N (/ = 0). 
The formal comoving resolution length is then L res - L/N - 
0.40/; _1 Mpc, corresponding to k res - 2nlL res =* 15.7 h Mpc" 1 . 
The lower panel of Figure Q]however shows an almost vanishing 
spectrum at this wavenumber. As a matter of fact, the spectrum 
is already smeared by a factor ~ 0.8 at k =* 3. Accordingly, quite 
in general, spectra are unreliable when k increases above 20 % of 
the formal resolution limit k res . In this case, clearly, / ought to 
be suitably increased. We then need to ensure that no shot noise 
will interfere with this. 

As a matter of fact, even a random distribution of N particles 
exhibits a non-vanishing spectrum oc N . Smith et al. (2003) 
performed a detailed analysis of this numerical noise, by study- 
ing the evolution of 236 3 particles initially distributed with a 
spectrum oc k~ 2 . 

Here, we define the scale where the numerical noise begins 
to affect spectra by evaluating a numerical derivative dP(k, z)/dk 
and fixing it to the k values where the spectrum turns from 
decreasing to increasing. This evaluation requires a suitable 
smoothing of the spectra (five-point Gaussian averaging) and 
the derivative itself includes an additional smearing (again of 
five-point Gaussian averaging) to appear as shown in Figure |2] 
Taking into account the different features of our simulations, our 
results are consistent with those of Smith et al. (2003). 

The behaviour of dP(k,z)/dk shown in Figure 2 then indi- 
cates a limitation of the use of spectra shifting from k =* 40- 
50 /zMpc~' (comoving L/A _1 Mpc - 0.09) at z = 0, to 10- 
12 (0.42-0.52) at z = 2.55 . As we show, the contributions 
to shear observables from higher redshift become quite small. 



Accordingly, at z — the numerical noise affects k values ap- 
proximately down to k € /5 . The P4(k,z = 0) spectrum is then 
substantially unaffected by numerical noise, which still causes 
some distorsion in P5, and limits the use of spectra at z - 0. 

A key issue that we anticipate is that shear spectra, resulting 
from an integration along z over variable k values, require the 
knowledge of P(k,z) at progressively smaller k as z increases. 
Since the integration involves passing from a spatial to an an- 
gular spectrum, the same angle & ~ \0(2n j €j\ subtends an in- 
creasingly larger scale A ~ [0(2n/k)] at greater distances. This 
effect occurs approximately in parallel with the downward shift 
in k caused by the impact of numerical noise. Thus, the spectral 
range we need to consider at higher z is mostly noise free, once 
it is so at z - 0. 

We however emphasize that, in recent works dealing with 
the effects of baryons on the power spectra, no n value above 
N was used to produce spectra. With a box of ~ 60 /r'Mpc and 
N - 5 12, as used by Rudd et al. (2008), we then have L res =* 0. 12 
and kres ~ 50. However, as shown in the lower panel of Figure 
Q] spectral information is scarce above k ~ 10 =* k res /5 , even at 
z — . This place severe limitations on our ability to evaluate 
Pij{() at large I . In our case, a similar limit is attained only at 
z ~ 2.5, while at lower z values, fluctuation spectra are reliable 
well beyond which k value. 



3.2. Dark matter simulation: Simulation spectra versus 
Halofit 

In Figure|3] we explore the spectral evolution across the redshift 
range z = 0^1 , by showing k 3/2 P(k,z). Power spectra for the 
DMo simulations (solid curves) are obtained with / = 5 and are 
also compared with HF predictions (dotted curves). 

In the top panel, we give an overall picture of all spectra used 
up to z ~ 4, in a restricted k range. In the bottom panel, we widen 
the k range from to 3, but restrict ourselves to five spectra. 

These plots illustrate that there is a considerable disagree- 
ment between the HF predictions and, those of the simulations, 
significantly exceeding the ±3 per cent approximation claimed 
in Smith et al. (2002fl 

At z ~ 3-4, in the range k ~ 1-2 h Mpc -1 , we find discrepan- 
cies exceeding 15-20 %. In general, we note that there is a sys- 
tematic lack of power in HF, with respect to simulations, even in 
the region of a non-linearity onset. Such a lack of power bursts 
above k =* 10/zMpc -1 . As we discuss below, this is however the 
regime where the baryon physics should be most important and 
the predictions of the DMo simulation less reliable. When AGN 
feedback is included, Van Daalen et al. (201 1) also noticed that 
baryon physics affects the matter distribution down to k ~ 0.3- 
IhMpc- 1 . 

A robust analysis should certainly be based on more model 
realizations. We however recall that a lack of power would be 
the natural consequence of a shortage of long-wave contribu- 
tions. The HF expressions, which have been calibrated using 
simulations within boxes of sizes never exceeding 240 /t~' Mpc 
(see Jenkins et al. 1998, Smith et al. 2003), could actually lack 
some power, leading to a delayed onset of non-linearity. 



1 Smith et al (2002), however, claimed to reproduce the results of 
their simulation with HF only in a restricted range of k and z, and not in 
general 
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Fig. 3. Spectral evolution, for DMo simulations, in the redshift 
range used to compile shear spectra. Spectra obtained with / = 5 
are shown. Simulations are also compared with the HF predic- 
tions. In the top panel, we show all spectra (multiplied by k 3 ^ 2 , 
to allow an overall comparison) in a restricted k range. Redshift 
values are shown along aside each curve. In the bottom panel, we 
select five redshift values and plot spectra over a wider k inter- 
val. Substantial discrepancies from HF are evident, even at low 
k values, at any redshift. 



3.3. Hydrodynamical simulations: total spectra 

We then consider the effects of baryon physics, starting from the 
analysis of the non-radiative (GH) simulation. 

In Figure|4] the total fluctuation spectra of the GH simulation 
are compared to DMo spectra. The pressure support of the bary- 



onic component partly inhibits the increase in the fluctuations 
in the total matter distribution, thereby reducing the amplitude 
of the power spectrum on scales below the corresponding Jeans 
length. In Figure|5] we present the same comparison for the CSF 
simulation, which includes the effect of radiative cooling and 
star formation. 

In the GH case, the reduction in the power spectrum am- 
plitude due to baryons gradually moves towards larger k values 
as one moves to lower redshift. The onset of white noise can 
obscure this effect, which is hardly visible above z ~ 1 . The 
average interparticle separation is different in the DMo and GH 
simulations, since two particle populations (one collisionless the 
other collisional), displaced by half grid cell, are used to sample 
density and velocity fields in the initial conditions of the hy- 
drodynamical simulations. Therefore, twice as small as number 
of particles in the DMo simulation justifies the correspondingly 
higher level of white noise reached above the Nyquist frequency. 
This effect must be kept in mind in any comparison between 
DMo and hydrodynamical simulations. 

Owing mostly to the larger range covered by the y-axis, it is 
even harder to appreciate here a slight lack of power in the CSF 
simulation at k ~ 1 /zMpc -1 . The spectra for the CSF simulation 
are then systematically above those of the DM-only case. This 
increase in the small-scale power in the CSF case is due to the 
contraction of halos induced by gas overcooling (see also Rudd 
et al. 2008), which alse removes the pressure support responsible 
for the inhibition of the fluctuaction growth in the GH case. 

Quite interestingly, in both GH and CSF simulations, we 
found the inverse effect around k ~ 1 h Mpc" 1 , where a slight in- 
crease in the spectral amplitude appears. This effect is more ap- 
parent in the upper panel of Figure [6] which provides a zoomed 
view. 

The slight excess power around k ~ 1 AMpc -1 in the DMo 
simulation has also been by analyses of the power spectrum in 
cosmological N-body and hydrodynamical simulations ( Van 
Daalen et al. 20110)- This feature occurs in a dynamical range 
where the effects of non-linearity start to appear. A similar effect 
was also outlined by Rudd et al. (2008), also using both a non- 
radiative and a radiative simulation, similar to our GH and CSF 
cases, but carried out with the adaptive mesh refinement (ART) 
code (Kravtsov et al. 1997). In addition Casarini et al. (2011b) 
found a similar effect in hydrodynamical simulations carried out 
with the SPH Gasoline code (Wadsley et al. 2004), using a box 
of 64 /z~'Mpc on a side. 

The size of the effect found here can be appreciated from 
Figure [6] While Rudd et al. (2008) found an effect as large as 
ten per cent, Van Daalen et al. (201 1) and Casarini et al. (201 lb) 
found a smaller effect, ranging from one to two percent, thus 
closer to the one obtained from our analysis. Furthermore, differ- 
ences also exist in the k range where the effect arises. It is quite 
difficult to draw a strong conclusion from this comparison about 
the analyses based on simulations covering different dynamic 
ranges, and using different implementations of radiative cooling 
and star formation. Nevertheless, we note that the largest feature 
of power inversion appears in the simulation carried out with a 
Eulerian grid-based code, with respect to the other SPH-based 
analyses. 
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Fig. 4. Power spectra of density fluctuations for the non- 
radiative (GH) hydrodynamical and DMo simulations. Results 
are shown for / = 5. Top panel: evolution between z — 3 and 
z — for the GH (solid curves) and DMo (dotted curves) sim- 
ulations. Redshift values are indicated at the side of each curve. 
Here, the small spectral differences are hardly visible, up to the 
onset of white noise at large-A;. Bottom panel: ratio of power 
spectra of the DMo to GH simulations. A slight excess of power 
in the GH simulation is visible up to £ =* 10/iMpc -1 . The in- 
version at higher k value is covered by the onset of white noise 
(see also Figure 0. Note the different levels of white noise in 
the two simulations, owing to the different numbers of particles 
used (see text). 



3.4. Hydrodynamical simulations: behaviour of the different 
components 

In Figure|7] we show the power spectra of the DM and gas com- 
ponents in the hydrodynamical simulations, normalized to the 
corresponding power spectrum computed for the DMo case. The 
three panels show results at z — 0, 0.58, 1.51 (from top to bot- 
tom). The upper and lower parts of each panel show results plot- 
ted using linear and logarithmic units on the y-axis, respectively. 

We outline a few points that are directly observable there. 
First of all in the GH case, the total spectrum is mostly above 
the corresponding spectrum of the DM component, while the 
gas spectrum is increasingly below that of the DMo simulation. 

In contrast, in the CSF case the total spectrum exceeds that of 
the DM, above k ~ 9-10 hMpc~ l . This is presumably due to the 
contribution of stars, which are highly concentrated within ha- 
los, thus providing a large amount of fluctuation power on small 
scales. The same feature, although appearing at slightly different 
values of k, was also found by Rudd et al. (2008), Van Daalen et 
al. (201 1), and Fedeli et al. (201 1). 

Furthermore, as we have already pointed out, hydrodynamics 
acts as a brake on the growth of overdensities on scales smaller 
than the Jeans length. Therefore, in the non-radiative GH case 
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Fig. 5. The same as in F|4]but for the comparison between the ra- 
diative (CSF) simulation and the DMo simulation. Top panel: the 
larger amplitude of the CSF power spectrum is clearly visible, 
well before the onset of white noise at large-£. Bottom panel: 
ratio of the DMo to CSF spectra. A slight lack of power in CSF, 
up to k a 10/tMpc -1 , is hardly appreciable. 
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2 In this paper we compare our CSF simulation with the REF one in 
Van Daalen etal 2011 



Fig. 6. A blown-up view of the lower panels of Figures |4] [5] to 
show the slight inversion of the ratio of the power spectra of 
DMo to hydrodynamical simulations, which takes place around 
k ~ 1 /iMpc~' ; an effect previously found by various authors (see 
text). The effect is smaller in CSF than in GH. 



the fluctuation amplitude of the gas component becomes smaller 
than that of the DM component. 

The presence of radiative cooling allows gas to dissipate en- 
ergy and sink to the centre of DM halos, thereby increasing 
their concentration by adiabatic contraction (e.g. Gnedin et al. 
2004, Sellwood & MacGough 2005). At z - 0, this causes an in- 
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Fig. 7. Spectral evolution of the different components in the GH 
and CSF simulations. For comments, see text. 



crease in the fluctuation power on small scales k >~ 40 /iMpc~', 
namely on a length scale ~ 150/r'kpc, close to the scale length 
where CDM and gas density profiles begin to diverge from each 
other. 

While, in Van Dalen et al. (2011), the total spectra were 
found to overcome CDM around the same k as here, the scale 
where the gas spectrum overcomes the DMo spectrum is greater. 

Secondly ,the evolution in redshift of the total power spectrum 
and its components allows us, first of all, to follow the progres- 
sive onset of the inversion in the ratio of DMo to either CSF or 
GH spectra, as described in section 3.3 . Logarithmic plots are 
appropiate way of following this process. 

At redshift as high as z = 1.51, when the £-range below 
10/iMpc~' is clearly distinguishable above numerical noise, we 
can discerne an inversion bump at k ~ 2, in the GH total spec- 
trum, whose amplitude attains ~ 0.5 %. At z = 0.58, the inver- 
sion interval has widened in GH, reaching an amplitude ~ 1 %, 
and is also visible in the CSF case. Altogether, however, the evo- 
lution of the GH/DMo ratio exhibits a waving development, with 
a power deficit at short wavelengths. In CSF, the inversion bump 
is instead isolated and does not exceed 1 %, even at z — . 

Our third remark is that there are discrepancies between our find- 
ings and Rudd et al (2008). These authors, found a magnification 
of the inversion in the presence of cooling and star formation. 
Furthermore, they were unable to distinguish the "waving" be- 
haviour in the GH case, in the counter-inversion interval when k 
is between ~ 2 and 15. On the other hand, our finding confirms 
Casarini et al. (2011b), Van Daalen et al. (2011), and Fedeli et 
al. (201 1) results. These discrepancies might arise from the dif- 
ferent codes used or from the too small box used in Rudd et al. 
(2008). 

Our fourth is that the high-£ behaviour at high z allows us to ap- 
preciate in detail how the different level of numerical noise arises 
in the different components. In particular, the amounts of power 
for CDM and GAS in either GH or CSF coincide with that for 
DMo, because each component share the size of the grid where 
initial conditions were set, even if particles are displaced by half 
a cell. In the hydro cases, the total spectrum noise exceeds the 
DMo one by the known factor. 

The fifth remark is that when cooling is allowed, we clearly note 
that gas loses energy more rapidly than CDM, and it then falls 
towards the centres of the forming halos. It is mostly because of 
gas accumulation on small scales that the total power increases 
at large k values. It is also interesting to note that this process 
scales down with redshift. This could partly be caused by the in- 
crease in the concentrations in halo profiles, which is also by the 
presence of stars. However, we also note that baryon physics has 
characteristic length and time scales that are approximately con- 
stant in physical space, while k is comoving and thus increases 
with decreasing redshift. 

We finally note that it is indeed significant that the results of 
Casarini et al. (2011b), in spite of the different resolution, box 
size, SNa cooling and numerical code, are consistent with those 
found here. 



3.5. Connecting linear and simulation spectra 

To construct fairly well-normalized shear spectra, density fluc- 
tuation spectra must extend down to k ~ 0.01 /zMpc -1 , or even 
less, in the full linear regime. For instance, some window func- 
tions Wj(z) imply that there is a significant contribution from 
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power spectra to shear spectra at z ~ 1-1.5. The typical wave 
numbers for which fluctuation spectra are known are then 



k^(l[T -T( Z )] 



Where t is the conformal time in the metric 



ds z 



a 2 (T)(d T 2 



drj 2 ) 



(2) 



(3) 7 



77 and a(f) being the elementary space interval and the scale 
factor, respectively, and t„ and t(z) are the conformal age of 
the Universe and its conformal age at redshift z . Accordingly, 
t - t(z) ~ 3000 h- l Mpc and I ~ 30 yields k ~ 10~ 2 /zMpc" 1 . 

In turn, this implies that using simulations within boxes of 
L ~ 60/!~'Mpc on a side (i.e. reaching k - 2n/L ~ 0.1 /iMpc -1 ), 
no direct connection between linear and simulation spectra can 
be made. In this case, one has then to resort to non-linear ap- 
proximations, such as HF, to cover an intermediate A:-range, up 
to a £ value where simulations provide a sufficiently large num- 
ber of Fourier modes (for L - 60/i~ 1 Mpc, discreteness effects 
in the sampling of these modes are significant even at k ~ 0.3- 
0.5/iMpc 1 ). Delicate normalization problems then often have 
to be solved. Although sharing the same linear o-$, DMo and HF 
spectra differ at most k values, as shown in Figure 3. Therefore, 
when a too small simulation box is used, the discrepancy often 
has an opposite sign, and exhibits a different size for different 
cosmologies. 

The alternative option of comparing simulations with equal 
values of non-linear erg improves the situation, since it over- 
comes at least the part of the discrepancy arising from the dif- 
ferent timing of growth. This however requires a complex "trial 
and error" procedure, which clearly becomes less and less de- 
manding as the box size is increased. For instance, Casarini et 
al. (2011a), by using L =* 250/i~ 1 Mpc, reach a satisfactory con- 
vergence with only a couple of iterations. 

Directly connecting the power spectra from simulations and 
linear theory eliminates any such problem. Although, at small 
k, simulation spectra exhibit a significant discreteness, our box 
size of 410/j _1 Mpc provides a connection between simulation 
and linear spectra in a A: range where a sparse sampling of k- 
modes induces only marginal effects. Following the same proce- 
dure as in Casarini et al. (2011a), we slightly smooth discrete- 
ness effects above k ^ 0.1 by averaging each P{k„) value with 
P(k„^i) and P(k n+1 ), giving ^_ lnn+1 P(^„-i,„,„ + i) (a selected so 
that kl^Pikn-y) and h% +1 P(k„ + i) are equal) weight 0.15, 0.7, 
0.15 , respectively. We verified that small variations of these 
weights have no visible impact on our final results. 

In Figure [8] we show the connection between the linear and 
the simulation spectra for the CSF simulation at z = 0. The linear 
spectrum shown here is the same as that used to generate our 
initial conditions of the simulation. 

We note here that the spectrum produced by the simula- 
tion still exhibits the expected BAOs at the correct k values. 
This is true even where discreteness effects are still dominant. 
A mild non-linearity is already visible above k ^ 0.1 /iMpc -1 . 
The BAOs expected up to k ~ 0.4/zMpc ~' can also however 
be traced in the simulation spectrum, despite the small amount 
of noise arising from sample variance. Such a small amount of 
contamination of the signal can easily be damped by running 
simulations of different realizations of the same model in a sim- 
ilar box and averaging among them. 




-1 1 

Log(k/hMpc) 

Fig. 8. Connecting linear and simulation spectra for the hydro- 
dynamical radiative (CSF) simulation at z - 0. The line colours 
are as follows: blue is the linear spectrum; dotted red is the sim- 
ulation spectrum; black is the connecting spectrum. The junction 
area (marked with the dottet rectangle) is magnified in the inset 
box, to also show the (mild) effect of the three-point smoothing 
on the simulation spectrum. Starting from large k, we pass from 
the simulation spectrum to the linear spectrum, where the former 
one approaches the latter for the first time. 

4. From fluctuation to shear spectra 

4.1. Theory 

The angular power spectra for the weak lensing convergence are 
given by the convolution of the matter power spectrum with win- 
dow functions that also account for the redshift distribution of 
the population of lensed galaxies. Following Hu (1999), if we 
bin the galaxies into n redshift intervals, the convergence power 
spectrum between the i and j tomographic beams, covering the 
redshift intervals A, and Aj, is given by 



PijM 



XT0 
duWi(u)Wj(u)P\-,u 



(4) 



Here Hq is the Hubble constant at present time, where h is in 
units of 100 kms~' Mpc" 1 , to is the conformal time correspond- 
ing to present cosmic age, and P(k, u) is the fluctuation spec- 
trum at the conformal time r — To — u . This angular conver- 
gence spectrum exhibits no explicit dependence on h, while it 
depends only indirectly on the cosmic expansion history through 
the Wi(u) window functions. In the literature (Hu 1999, Casarini 
et al. 201 la), n - 1, 3, or 5 bins have been considered. Here we 
mostly refer to a five-bin case and the bin limits z, are selected 
so as to ensure that there is the same number of galaxies per bin. 
We then define the window functions W,(m), which tell us 
how clustered matter acts on the galaxies in the j-th bin. We 
assume that the number density of galaxies per unit redshift in- 
terval and solid angle is given by 



n(z) 



dQ.dz \Zo 



exp 



(5) 
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Fig. 9. Upper panel: Dependence on u of the P(£lu, u) integrand 
factor, for I = 200, 400, 1000, ,3800. The curves are ob- 
tained by plotting the 10,000 points used in Riemann integration. 
Lower panel: window functions W,(m) in the five-bin case (u in 
Mpc). 



with 
C = 



B 



*rffl] 



(6) 



Here we take the usual values of A - 2 and B = 1.5, so that 
C = 1.5/zo where z„, - 0.9 and zo - z m / 1.412 (z,„ median 
redshift; see, e.g., Refregier et al. 2006, La Vacca & Colombo 
2008 and references therein). 

This distribution is then considered within the limits of the 
redshift bins. Following Hu et al. (2006), we account for the dis- 
crepancies between the photometric redshifts, on which the red- 
shift distribution is based, and the true redshifts, by defining the 
filters 



n,(z) : 
_ 1 

~ 2 



f 



Erf 



dz' 



1 



V2^cr( Z ) 



exp 



iz-z') 2 
2cr 2 (z) 



ZphJ 



V2o-(z) 



I - \ ^ A Z ph,i 



-Erf 



\ V2cr(z)/_ 



(7) 



where cr(z) = 0.05 (1 + z) (see, e.g., Amara & Refregier 2007, 
for the motivation of this parameter choice). We point out that 
this expression for the filters assumes a Gaussian distribution for 
the errors in photo-z's. An accurate calibration of this error dis- 
tribution is clearly of paramount importance to avoid introduc- 
ing biases into the reconstruction of power spectra derived from 
weak lensing tomography (e.g. see Bernstein & Ma 2008). 
If we define 



Di(z) = n(z)n,(z) , 

we can then introduce the distributions 
A(z) 



Si(z) 



f^D^W 



(8) 



(9) 
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Fig. 10. Assumed redshift distribution of galaxies and their ex- 
pected partition among the five bins. In the inner panel we plot 
u (in Mpc) versus z. 

which indicate the actual redshift distribution of the galaxies be- 
longing to each bin (see Figure ITOb. From them, we derive the 
filter functions 



Wi(z) = ^mFKzXl+z), 



(10) 
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Log[k=*/(T -T)] 



Fig. 11. Domains of integration of the matter power spectrum 
P(k, t) on the k,r - t plane, to compute the convergence spec- 
trum. Each solid line corresponds to a value of I ranging from 
I = 4000 to / = 40000 with increments of M = 4000. The two 
horizontal dashed lines mark the value of the conformal time 
corresponding to redshifts z = 0.1 and 0.2. 
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Fig. 12. Numerical noise limit (magenta band) compared with 
the maximum k needed to calculate the shear spectra at various 
I values (aside the curves). The k values above which shot noise 
exceeds the physical spectrum is approximated here by the ex- 
pression log(fc//;Mpc _1 ) ~ 1.6 — 1.05 log(l + z), in accordance 
with the results of the method described in Sec. 3. 
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Fig. 13. Dependence on u — tq - r of the integrand in eq. |@). 
The case i — j — 1, which has the largest contribution at low 
redshift, is shown. The function is shown for I — 4000 (lower 
panel) and 40000 (upper panel). We note the different ordinate 
scales and the 1ow-m cut in the latter case, where we estimate a 
reduction of P u (40000) by ~ 4 per cent. 

where information about the galaxy redshift distribution is con- 
tained in the quantity 



F;(z)= f dz' 6i{z')u{z,z!)lu{z') . 

Jhzi 



(11) 



In the above expression, u(z, z') is the (positive) conformal time 
distance between z and z' ■ In Figure [9] we show the W, profiles 
in the five-bin case; they are used in eq. ©, by replacing the 
dependence on z with the dependence on the conformal time- 
distance u — to - t (see the inner box in the very Figure [9}- 
Figure |9]highlights that in order to predict shear spectra we need 
fluctuation spectra P{k,z) up to z ~ 2.5, i.e. u =* 5000 Mpc. In 
turn, this means that tomographic shear spectra for the assumed 
redshift distribution will hardly provide sufficient information 
about the evolution of the matter density power spectrum at z >~ 
2 . The redshift range z =0-2, however, is the one where one 
expects to be able to see the effect of the presence of DE. 

The range of k values over which we need to sample the 
matter power spectrum is illustrated in Figure [TT] This figure 
shows how we can derive shear spectra up to t a* 40000 (and 
beyond), where we expect to have underestimated them by $ 4 
per cent, as we argue below. 

The point is not to run the risk of using regions of the fluc- 
tuation spectrum where numerical noise is significant. In Figure 
fT2l we then compare the numerical noise limit (magenta band) 
estimated from power spectrum with f = 5, with k = lju, for 
I = 2000, 10000 and 50000. The figure confirms that a danger 
might exist at small u values. However, even for t - 50, 000 the 
risks concern spectra at z <~ 0.13 and can be directly avoided 
by using spectra with / = 4, which, at such low z, are cut- 
off before the onset of numerical noise. Rather than an over- 
estimate, we then have an underestimate of Ptj(€) for large I 
values. The underestimate is however small. For small u val- 
ues, the functions Wi(u)Wj(u)P(l/u, u) provide negligible contri- 
butions, as the window functions remains essentially constant, 
while the spectrum decreases by several orders of magnitude, 
towards large k's, yielding a rapid decrease with u. This in- 
dependence is also shown in the upper panel of Figure [9] for a 
sample of small t values. 

In Figure [T3] we then show the overall effect of the high- 
k spectral cutoff for i,j = 1, the spectrum most sensi- 
tive to low-z contributions, by plotting the whole integrand 
W\(u)W\(u)P{llu,u) for I = 4000 and 40000 . For i = 4000, 
the function appears quite regular, while for C = 40000, we have 
a 1ow-m cut-off, arising from the spectral cut-off. This cut leads 
to an underestimate of Pi\(l), starting at t ~ 25000, reaching 
~ 4 % at t - 40000 (the case plotted) and attaining ~ 10 % for 
t - 50000. Deficits are smaller at larger i, j values. 

4.2. Tomographic shear spectra 

Before using eq. © to evaluate the shear spectra up to large I 
value, we proceed with a few further tests. 

In Figure [14] we show the effect of artificially cutting 
the spectra at a redshift-independent value k cut . When k cut - 
50 /iMpc -1 , the logarithmic plot of Pi i does not appreciably dif- 
fer from the one computed using full spectra. This plot can be 
compared with Figure Q3] which shows the effects of different 
resolutions in the fluctuation spectra. This confirms that lower 
resolution grids, yielding spectra that decline systematically well 
before k ~ 40 /zMpc -1 , are inadequate to reach large I values. 

In Figure[l6] we then show the shear spectra P,,(0 and Pu{() 
for i = 1, ..., 5, for both the hydrodynamical CSF simulation and 
for DMo iV-body simulation. The extra power at large k's in 
the spectrum of the CSF simulation is clearly visible at large I 
(> 1000). The relative difference in the power spectra between 
these two simulations can be more clearly appreciated in Figure 
[P71 This figure exhibits an interesting feature: at t ~ 10 2 -10 3 , 
there is a short interval where the shear spectra of the CSF sim- 
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Fig. 14. Effects on P\ \{E) of cutting the fluctuation spectra for the 

CSF simulations at different values of k cut . Progressively lower Fig. 16. Tomographic shear spectra in the five-bin case. In the 

curves correspond to increasing values of k cut , as reported by the upper (lower) panel, we plot P,,(7) (Pi,(/)) with i - 1 , ..., 5. for 

labels. the hydrodynamical CSF and DMo N-body simulations. 
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Fig. 15. Large-^ cutoff of Pi i {(), obtained by integrating spec- 
tra obtained with different grids. The effective number n of grid 
points used is indicated in the frame. Note that the spectrum for 
the DMo simulation computed for n - 16384 is lower than the 
spectrum of the CSF simulation computed for n - 8192. 



ulation are slightly below the corresponding ones of the DMo 
simulation. This interval is more evident in Pjj(t) for low i, j, 
and is a consequence of the power inversion in fluctuation spec- 
tra (see Figure 6). We point out however that the difference be- 
tween shear spectra is dominated by the large shift from P(k) 
at k >~ 2-4 . When integrating over u, such shifts affect all I. 
Hence, the detailed features of Figure [6] are diminished by the 
large-fe shift. When comparing GH and DMo, these features ap- 
pear more clearly (see below). 

On the basis of the results shown in Figure 3, one may won- 
der to what accuracy HF can reproduce simulation results for 
the shear spectra. This is shown in Figure [18] which highlights a 
major difference, that is sizable even at t - 1000. This confirms 



Fig. 17. Relative difference in the shear spectra for the CSF and 
the DMo simulations, for the three different values of i, j re- 
ported in the labels. The sign of the shift changes in two points. 
Between them, the logarithm of the modulus is plotted. 



that using HF to deal with non-linear shear predictions can lead 
to substantial inaccuracies. 

We finally test the effects of considering the effect of non- 
radiative hydrodynamics on the shear spectra. In Figure [19] we 
compare Pu(€) spectra for the GH hydrodynamics case with N- 
body spectra, as done for CSF spectra in the top frame of Figure 
[T6l We find that any difference between DMo and GH spectra is 
small, never exceeding 10 % up to £ ^ 50000. An interval where 
DMo spectra (marginally) are brighter than those of GH exists, 
between I - 1000 and ~ 10000. This inversion ceases at a higher 
( value for greater values of the indeces i. Altogether, this plots 
reproduces, in terms of {, features that can be discerned when 
comparing P(k) spectra, in Figure [6] 

Before concluding this section, we comment on the accu- 
racy of the numerical integration yielding the shear spectra. 
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Fig. 18. N-body shear spectra compared with shear spectra ob- 
tained using HF, for the same model. 




Fig. 19. Shear spectra from GH and DMo simulations nearly 
overlap. In the lower frame, we then plot the shift (Pu(ptfo) ~ 

Pii(GH))/Pn(CH), which becomes significant only above I =* 1000, 
as expected, and but remains within 10%, being smaller for 
larger i values. 



Numerical inaccuracy could cause a dangerous misunderstand- 
ing about the transfer of fluctuation spectra features onto shear 
spectra. As a matter of fact, all Pjj(l) shown so far were ob- 
tained by performing the integration in eq. (2) using a simple 
Riemann algorithm, by summing over N p ^ 10 4 integration 
points evenly distributed between u — and 6500 Mpc (z - 
and 3). For low i, j, quite a few points are doomed to fall into 
a region where the window functions ensure that the integrand 
is negligible. For instance, Figure [9] shows that W\ ~ 10~ 4 at 



z ~ 0.5 (u ~ 1500 Mpc); even in this worst case, about ~ 2400 
points are then significant. Convergence was however numeri- 
cally tested, by considering a variable N p , and is already satis- 
factory for N p = 2000 . 



5. Conclusions 

We have presented results on the computation of density fluctu- 
ation power spectra and tomographic cosmic-shear power spec- 
tra, from a set of A^-body and SPH hydrodynamical simulations 
of a concordance ACDM model. Simulations were carried out 
in a box with a comoving size of 410/i~'Mpc using N = 1024 3 
dark matter (DM) particles and an equal number of gas particles 
in the hydrodynamical simulations, with a Plummer-equivalent 
force resolution of 7.5 ft kpc. This enabled us to compute spec- 
tra up to k - 500/zMpc , whose reliability is only limited by 
numerical noise effects that, even at z = 0, are still evident for 
k> (50-70) ftMpc -1 . 

In addition to a DM-only N-body simulation, two hydro- 
dynamical simulations were carried out, the first one based on 
simple non-radiative physics, and the second one including ra- 
diative cooling, star formation, and the effect of SN feedback in 
the form of galactic winds. Since hydrodynamical simulations 
are based on twice as many particles, their power spectra be- 
ing characterized by a relatively low level of high-A: white noise, 
with respect to the DM-only simulation. 

Taking advantage of the fairly large size of the simulation 
box, we performed a direct connection between spectra from 
simulations and linear-theory, in a £-range where the latter is ex- 
pected to hold to good accuracy. A comparison with the predic- 
tions of HF confirms the unreliability of this algorithm in provid- 
ing precise predictions for non-linear power spectra. The deficit 
of power from HF appears to be more and more severe at higher 
redshifts, likely due to the delayed onset of non-linearity pre- 
dicted by this approach. When non-linearity is finally fully de- 
veloped, discrepancies becomes smaller, although are still quite 
significant. This discrepancy could possibly arise from the lack 
of long-wave contributions to the spectra used to build HF, that 
were obtained from simulations in boxes whose side never ex- 
ceeded ~ 240 ft^'Mpc. One cannot however exclude that the re- 
alization of this simulation is notably above average. 

This finding highlights the need to study the non-linearity 
onset by using large enough boxes. In a number of recent papers 
studying the effect of baryon physics on the non-linear power 
spectrum (e.g., Rudd et al. 2008, Van Daalen et al. 201 1), boxes 
with sides of mostly <9(60ft~'Mpc) were used. Van Daalen et 
al. (2011), also used a slightly larger box of 100ft~'Mpc on a 
side. In our opinion, using boxes of such limited size can lead to 
non-numerically converged estimates of the non-linear power 
spectrum. 

In this paper, we have not debated the nature of DE, and 
referred to purely ACDM cosmologies. Moreover, AGN feed- 
back was not included, at variance from both Van Daalen et al. 
(2011) and Semboloni et al. (2011), who however made use of 
quite a smaller box. Because of spectral discreteness, to obtain 
shear spectra, they ought then to make recourse to approximated 
spectral expressions, even before the spectral scale approaches 
the box size, being there still far from a fully linear regime. The 
whole spectral normalization may then become imprecise, par- 
ticulary when spectral expressions suitable to ACDM models are 
tentatively extended to other cosmologies (see, e.g., Casarini el 
al. 201 la). In contrast, ournumerical spectra match a scale range 
where linearity is substantially unaffected (see Figure [8). The 
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box size is therefore a critical issue, as no approximated expres- 
sion is needed in our analysis. 

As a matter of fact, following the onset of non-linearity with- 
out introducing any source of bias, may become a key issue when 
tomographic shear spectra of different DE models are compared. 
The models often differ mostly in terms of the timing of the non- 
linearity onset. When lacking the correct normalization between 
spectra from simulations and linear theory, this timing is dis- 
torted and model comparisons can be biased. 

As in our study the effect of introducing baryons, our anal- 
ysis confirms that significant differences from A^-body results 
appear at large k values when radiative physics is included. In 
this case, the increase in small-scale power originates from the 
sinking of a significant amount of cooled baryons in the central 
regions of DM halos (e.g., Gnedin et al. 2004, Jing et al. 2006, 
Rudd et al. 2008, Duffy et al. 2010, van Daalen et al. 201 1). 

As for the computation of a tomographic shear spectrum, we 
highlight that the large size of the box size used has allowed us to 
obtain them from the density fluctuation spectra by directly re- 
lating the linear spectrum to that measured from the simulations. 
In most previous analyses, which are based on smaller simu- 
lation boxes, HF expressions were used to cover a fairly wide 
intermediate-range of k values. The spectra Pij{€) could then 
be computed for I values reaching 50000. The simulations used 
here allowed us to achieve a precision better than 1 per cent for 
I <~ 25000, with an underestimate of the spectrum amplitude of 
within 10 per cent at t ~ 50000, in the worst case, thanks also to 
pushing numerical noise interference up to k > 40-50 /iMpc" 1 . 

Our results highlight the importance of accurately calibrating 
the subtle effects of the propagation of signal from fluctuation 
spectra to shear spectra, and the relevance of detailed large-scale 
simulations in tracing the transition from linear to non-linear 
scales. 

ACKNOWLEDGMENTS. LC acknowledges the Brazilian re- 
search Institutions FAPES and CNPq, and the Observatory 
of Paris-Meudon for their financial support. SAB acknowl- 
edges the support of CIFS. We acknowledge partial support by 
the European Commissions FP7 Marie Curie Initial Training 
Network CosmoComp (PITN-GA-2009-238356), by the PRIN- 
INAF09 project Towards an Italian Network for Computational 
Cosmology, by the PD51 INFN grant. We are grateful to 
Volker Springel for making the non-public GADGET-3 code 
available to us. Simulations were carried out at the CINECA 
Supercomputing Centre in Bologna, with CPU time allocated 
through a ISCRA proposal. 



Heitmann K., White M., Wagner C, Habib S., Higdon D., 2010, ApJ, 715, 104 
Hilbert S., Hartlap J., White S.D.M., Schneider P., 2009, A&A, 499, 31 
Hoekstra H., Mellier Y., Van Waerbede L., Semboloni E., Fu L., Hudson M.J., 

Parker L., Tereno I., Benabed K., 2006, ApJ 647, 1 16 
HuW., ApJ522L21 (1999) 
Hu W., Ma Z. and Huterer D., ApJ 636, 21 (2006) 
Huterer D. and Takada M., 2005, Astropart. Phys. 23, 369 
Klypin A. & Holtzman J., 1997, astro-ph/9712217 
Kravtsov, A. V., Klypin, A. A., & Khokhlov, A. M. 1997, ApJS, 111, 73 
Jenkins A., Frenk C.S., Pearce F.R. et al , 1998, ApJ 499, 20 
Jing Y.P., Zhang P., Lin W.P., Gao L., Springel V, 2006, ApJ 649, LI 19 
Komatsu E., Smith, K. M., Dunkley, J., et al. 201 1, ApJS, 192,18 
Laureijs, R., Amiaux, J., Arduini, S., et al. 201 1, arXi v:1110.3"T93l 
La Vacca G. & Colombo L., 2008, JCAP 0803.1640 
Massey, R., Rhodes, J., Leauthaud, A., et al. 2007, ApJS, 172, 239 
Refregier A. et al. 2006, Proceeding of SPIE - Astronomical Telescopes & 

Instrumentation, Orlando, May 2006, & astroph/06 10062 
Refregier A., Amara A., Kitching T. D. et al. 2010, |arXiv:1001.0061 1 
Rudd D.H., Zentner A.R., Kravtsov A.V., 2008, ApJ 672, 19 
Schrabback T. et al., 2010, A&A 516, 63 
Sellwood, J. A., & McGaugh, S. S. 2005, ApJ, 634, 70 
Semboloni E., Hoekstra H, Schaye J., van Daalen M.P., McCarthy I.G, 2011, 

arXiv: 1105. 1075 MNRAS (in press) 
Seo H.-J., Sato M., Takada M., Dodelson S., arXivrl 109.5639 
Smith R. E., Peacock J. A., Jenkins A. et al. 2003, MNRAS 341, 131 IS 
Springel, V. 2005, MNRAS, 364, 1 105 
Springel, V, & Hernquist, L. 2003, MNRAS, 339, 289 
Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253 

Tornatore, L., Borgani, S., Dolag, K., & Matteucci, F. 2007, MNRAS, 382, 1050 
Van Daalen M.P., Schaye J., Booth CM., Dalla Vecchia C, 2011, 

larXiv:1104.lT74 MNRAS (in press) 

Viel M., Markovic K., Baldi M., Weller J., 2011, arXiv: 1107.4094 
Wadsley, J. W., Stadel, J., & Quinn, T. 2004, New A, 9, 137 

Appendix A. Computation of power spectrum with folding 

We consider a one-dimensional distribution of N particles on a 
segment of length L, i.e. with abscissas < x r < L , given by 



A(x) = y S(x - x r ) ■ 
Its Fourier transform is 



r=\ 



P(k„) 



1 r L N 1 N 

VI Jo tf vT tf 



(Al) 



1ninx T jL 



(AT) 
with k„ = y' ■ «■ We then evaluate P(k) using a simple FFT 
algorithm, obtaining the spectrum up to « = 2^ (fi positive in- 
teger). For instance, we consider fi - 9, yielding h — 512 
and therefore allowing a formal resolution down to the scale 
/l = L/512 = L/2 9 . 

From A(x), we can obtain another distribution 
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A2(x) = Y^S( X -4)= Y, $(X-Xr)+ J] $X~(Xr-L/2)], 

r=\ x r <L/2 x r >L/2 

(A3) 
whose particles are all set on points x' r < L/2, either because they 
are originally so, or because their abscissa has been lowered by 
L/2 . 

Its Fourier transform is given by 



P(k m ) = 



1 



\ 1 2n i ImXyjL 
L.v,<L/2 



z 



+ 2, ' 

x r >L/2 



2ni2m(x r -L/2)/L 



(A4) 



with k m - jj2 ' m ■ The second term in the r.h.s. however also 



reads 

V^ 2ni2m[x r -L/2]/L 



Xr>L/2 



-z- 

x r >L/2 



27ii2mx r /L~2n i 



im _ \ 1 2n i2mx r /L 
x r >L/2 
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as exp[27r im] = 1 for any integer m. 

Accordingly, the spectra of A and A2 coincide, for n = 2m, 
apart of a factor a/2. We again choose m — 2^ with jj. — 9. The 
formal resolution scale is then (L/2)/512 = L/2 10 . Accordingly, 
we double the resolution, although we evaluate the spectral har- 
monics only for even n values. 

The generalizations of this procedure to L/4, L/8, etc., and 
then three dimensions are straightforward. 
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